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I. INTRODUCTION 

During the past decade several new charmonium reso- 
nances were discovered, primarily by experiments at the 
two B-factories but also by CLEO-c and at the Teva- 
tron. With BES-III and the LHC collecting data, pos- 
sible Super- i? factories and the planned PANDA experi- 
ment [l| at the FAIR facility, experimental prospects to 
study these states in more detail and to discover further 
resonances are very promising. For an overview, see e.g. 
Refs. 

Current phenomenological debates focus on the X, Y 
and Z resonances that are close to or above open charm 
thresholds. At least four different frameworks have been 
suggested to accommodate these states: 

• D'-*^D molecules (or deusons) [8l-[l3l|. composed of 
a charmed meson and antimeson D, 

• tetraquark states (or baryonia) [l^ - |T8| consisting of 
diquark-antidiquark pairs, bound by QCD forces. 



• ccg hybrid (or hermaphrodite) states |19l - l22| con- 
sisting of a charm-anticharm quark pair and addi- 
tional gluons, and 

• a compact cc core, bound inside a light meson, 
hadro-charmonium (23l . [2^ . 

One example of a molecule or tetraquark candidate is the 
A'(3872). The F(4260) can at present be accommodated 
as a hybrid or as a hadro-charmonium state while the 
Z+(4430) (if confirmed) could either be a molecule or 
hadro-charmonium. 

The new states also pose novel challenges to lattice 
simulations. In the case of standard charmonia that can 
be classified according to a non-relativistic quark model, 
the sizeable quark mass rric > A, where A denotes a typ- 
ical hadronic binding energy, represents the main diffi- 
culty: lattice artefacts, that in our case are of OKmca)"^], 
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are usually not small at currently available lattice spac- 
ings a. In the T case the b quark mass can be integrated 
out and an effective field theory, non-relativistic QCD 
(NRQCD), simulated on the lattice [2^, H^. However, 
the charm quark mass rric is not sufficiently large to al- 
low for this. In this case higher order perturbative or 
non-perturbative corrections will be sizeable. Therefore 
the charm quark needs to be simulated using a relativistic 
action. 

One would expect observables that are very sensitive 
to the mass tUc to be more strongly affected by lattice 
artefacts than those that are insensitive to the precise 
value of this mass. Using effective field theory methods 
like the Fermilab approach to heavy quarks l27H29i or 
NRQCD [30^ and potential NRQCD (pNRQCD) [Si HI 
some insight can be gained into this. For instance, charm 
quark mass effects on spin-averaged splittings are sup- 
pressed by a factor of the squared average relative veloc- 
ity of the charm quarks v'^. Momentum exchanges oc rricV 
in turn become relevant for the finestructure. Finally, 
lattice spacing effects on determinations of the mass pa- 
rameter rric from the charmonium spectrum are not sup- 
pressed by any powers of v. This means that a computa- 
tion of the spin-averaged spectrum will be less demanding 
with respect to the continuum limit extrapolation than 
predictions of the finestructure or of the charm quark 
mass. 

The standard spectroscopy of charmonium states in- 
cluding the continuum limit extrapolation is well un- 
derstood in the quenched approximation to QCD, see 
e.g. [33I and references therein, and several new re- 
sults including sea quark flavours exist, on isotropic lat- 
tices 29, 34 37] as well as on anisotropic lattices [l^ 
that employ a smaller temporal than spatial resolution 
at ^ fls, to lessen the severity of the scale separation 
rricV > rricv'^. 

An accurate reproduction of the charmonium finestruc- 
ture in the continuum limit represents an important test 
of QCD and of lattice methods. However, taking the con- 
tinuum limit may be less vital to reproduce qualitative 
features of loosely bound open charm threshold states 
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that are spatially more extended. In this case one needs 
to consider the mixing of states created by two-quark 
and by four-quark operators. Some pioneering studies 
have already been done, creating states with cqqc op- 
erators ^33, i40j where q denotes a light quark flavour. 
However, so far disconnected quark loop diagrams and 
hence annihilation channels have been neglected. More- 
over, lattice studies of the light quark sector, see e.g . 
Ref. [4l[, and of string breaking in the static limit [l^l 
teach us that these diagrams and, in particular, mixing 
between states created by cc and cqqc operators can be 
important. 

Here we will explore methods needed to systematically 
study charmonium threshold states and apply these to 
phenomenology. This article is organized as follows. In 
Sec. [H] we will introduce our methods, namely the gauge 
ensembles that are being used, the smeared operators 
that enter our variational analysis and the all-to-all prop- 
agator techniques. In Sec. Illll we present results on stan- 
dard charmonium spectroscopy at a finite lattice spacing, 
employing two-quark (cc) creation operators only. This 
includes higher spin and exotic states and provides us 
with the improved operators that are needed in Sees. IIVI 
andlVl In Sec. IIVI we investigate the mixing between cc 
and qq operators. This will yield an upper limit to the 
mixing between the flavour-singlet rjc and 77 mesons that 
in principle could have an effect on the S-wave finestruc- 
ture. Finally, in Sec.|V]we investigate the contribution of 
four-quark ccqq components to radially excited charmo- 
nium states, before we conclude in Sec. lVIl Some prelimi- 
nary results were presented at lattice conferences [4a - |45| . 



II. SIMULATION DETAILS AND METHODS 



TABLE I. Identifier (ID), simulation parameters, charm 
quark K-value (Kcharm) and the number of analysed effectively 
statistically independent gauge configurations of our runs. 



ID 


P 


K 


volume 


mps/GcV 


a/fm 


L/fm 


>i^charm 




® 


5.20 


0.13420 


16^ 


X 32 


1.007(2) 


0.1145 


1.83 


0.1163 


100 


® 


5.29 


0.13620 


24^ 


X 48 


0.400(1) 


0.0770 


1.84 


0.1245 


130 


(3) 


5.29 


0.13632 


24^ 


X 48 


0.280(1) 


0.0767 


1.84 


0.1244 


100 



physical light pseudoscalar mas^ [13], tops = in^^^. 
The measured values of rg /a not only depend on the in- 
verse lattice coupling /3 but also on the mass parameter 
K. One can now decide to define a lattice spacing a(/3, k) 
or replace this by a chirally extrapolated a(/3). After per- 
forming a chiral extrapolation in the sea quark mass, the 
results of the two choices obviously should agree for phys- 
ical observables. Since in this exploratory study we do 
not attempt such an extrapolation, we decide to set the 
lattice spacing from the ro/a(/3, k) values, as determined 
at the investigated sea quark k parameters. 

This leaves us with the charm quark mass as the only 
free parameter which we set by tuning, 

'^15 = \ i^va + 3toj/v[,) , (1) 

to its experimental value of jsj, (3067.8 ± 0.4) MeV. 

The ensembles and (2) are used to optimise the 
smearing functions. Our study of mixing between the 
rjc and the light quark ?7-meson is performed on where 
the mass gap between these states is smallest so that 
one may expect the biggest effect. For the mixing with 
threshold states ensemble (J) is used because in this case 
light D-meson masses are mandatory. 



A. Gauge configurations 

We base our study on np — 2 configurations of 
the QCDSF Collaboration generated using the non- 
perturbatively improved Sheikholeslami-Wohlert (clover) 
Fermion action [46] and the Wilson gauge action, pro- 
vided b y th e QCDSF collaboration. Details can be found 
in Ref. [4 71 ] . The charm quark mass to^ is not sufRciently 
heavy to allow for a non-relativistic action with control- 
lable systematics. Therefore, we use the same action for 
the charm quark as for the light sea/valence quarks, with 
a well-defined 0{a) improved continuum limit. Note that 
except for the value of the coefficient csw = cb ~ ce the 
clover action is identical to the version of the Fermilab 
heavy quark action used, e.g., in Ref. [2^ and our results 
at a finite lattice spacing a may be interpreted accord- 
ingly. 

We list the ensembles that we employ in Table HI to- 
gether with an identifier. The lattice spacing is set from 
the value rg ~ 0.467 fm. With this choice the nucleon 
mass agrees with experiment when extrapolated to the 



B. Tlie variational metliod 

We extract energy levels En from the decay of two- 
point Green functions in Euclidean time, 

a,{t) = mt)o]io)) (2) 

= ^«>;*e-^"*, (3) 

n>l 



A recent re-analysis yields somewhat different ro/a- values l4Sl . 
in particular at small quark masses. Here we ignore these de- 
velopments. Otherwise we would have to rerun all simulations, 
re-adjusting the charm quark mass by —5%, —6% and +8%, 
on ensembles (T), (2) and (J), respectively. However, most of the 
charmonium mass is given by 2mca so that only mass splittings 
will be affected by such a re-adjustment. Fortunately, the spin- 
averaged splittings were found to be rather insensitive to the 
charm quark mass 1491 while the main systematics regarding the 
finestructurc are the unrealistic sea quark content and the miss- 
ing continuum limit extrapolation. 
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TABLE II. Continuum spins that contribute to a given lattice 
representation. 



irrep. dimension continuum J 

Ai 1 0,4,- •• 

As 1 3,. . . 

E 2 2,4,... 

Ti 3 1,3,4,. . . 

Ta 3 2,3,4,. . . 



TABLE III. The "inverse" of Table [TTl Lattice spins that are 
"embedded" within each continuum spin. 

J Oh rop. dimensions 

1 

1 Ti 3 

2 E, T2 2+3 

3 ^2, Ti, Ta 1+3+3 

4 Au E,Ti, Ti 1+2+3+3 



where = {O\0i\n). In the case of the clover action 



link reflection positivity is violated and so in principle 
the above spectral decomposition vifith positive real en- 
ergy eigenvalues only becomes valid for sufficiently large 
Euclidean times. In practice for i > a we do not de- 
tect any violations. Oj are operators creating states of 
an isospin /, charm numbeiQ C, a given momentum and 
S0(3) ® Z2 ((8) Z2) J^^*^) quantum numbers. Note that 
on the lattice the infinite dimensional 0(3) group is bro- 
ken down to its octahedral Oh subgroup of order 48, with 
only ten {Ai, A2, E,Ti and T2 times parity) irreducible 
Bosonic representations. The mapping between the con- 
tinuum J spins and these Oh spins is given in Tables HIl 

-m 

The expectation value Eq. ([2]) will depend on the time 
difference between creation and destruction of the state 
so that for convenience we have set the source time to 
zero. Obviously, Cij{t) = C*^{t) is Hermitian and in 
our case we will use operators with phases so that it is 
real and positive definite ior t > a. For large times t 
the exponential decay of the Cij{t) entries will be gov- 
erned by the ground state energy i^i, or, for a momentum 
p = 0, by the ground state mass. Due to the transla- 
tional invariance of the expectation value, it is sufficient 
to perform this momentum projection at the sink. We do 
this for the standard spectroscopy so that we only need 
to generate point-to-all propagators in this case. Note 
that we still have the symmetry Cij(t) — Cji{t) in the 
limit of infinite statistics, however, the statistical errors 
of Cij{t) and Cji{t) for i ^ j will not be of similar sizes. 



^ Here we do not consider strangeness, beauty etc.. 



Replacing off-diagonal elements so that more smearing 
iterations (see Sec. IIIC|) are applied at the source than 
at the (momentum-projected) sink reduces the statistical 
errors. 

The convergence in Euclidean time of effective masses, 

C it) 

m,jMt + «/2) = In ^J^^^^ , (4) 

towards the ground state mass is affected by the quality 
of the ground state overlap q — \v}\'^ = |(l|6||0)p of the 
operator Oi. Having many different such operators at 
our disposal enables us not only to determine the ground 
state energy at small t-values where statistical errors are 
small but also allow us to access excited states, using the 
variational approach (sH - Iss! ]. also known as the general- 
ized eigenvalue approach. 

We choose a basis of operators Oi, i = 1, . . . ,N, de- 
stroying a colour singlet state within a given lattice rep- 
resentation. These operators may differ for example by 
their spatial extents or their Fock structures and they are 
usually not mutually orthogonal. These are then used to 
construct the correlation matrix Eq. ([2]). We now solve 
the symmetrized eigenvalue problem, 

C-'/^{to)Cit)C-'/^to)r{t,to) = X^itMWitM)- 

(5) 

Note that C-'^/^{to)C{t)C-^l'^{to) = 1 at the normal- 
ization time t — Iq: everything is expressed relative to 
the eigenbasis of C(to)- We order A^(t) > \^{t) > ■ ■ ■ > 
A^(t) > at large t. To ensure consistency over jack- 
knife samples, in the statistical analysis we also monitor 
the directions of the eigenvectors. Note that the orig- 
inal non-symmetrized definition of Ref. [Flj yields the 
same eigenvalues but different, non-orthogonal eigenvec- 
tors, 0"(i,to) -C-5(to)^"(t,io), 

C-\to)Cit)^"it, to) = A"(i, to)rit, to) . (6) 

If we choose to overly large then excited states will 
have died out in Euclidean time and the rank of C{to) 
will not be maximal, within the given statistical errors. 
For to chosen too small, C{t) will receive contributions 
from more than the N lowest lying states, resulting in 
unstable eigenvectors and eigenvalues. It can be shown 
that the eigenvalues behave like [53l |. 

A"(i, to) oc e-(*-*°) ^" [1 + 0(e-(*-*«) ^^-)] , (7) 

where Ai?„ is the energy difference between the energy of 
the first state not contained in the operator basi^, -Bat+i 



^ At least to first order in perturbation theory. To second order 
states with energies < En can contribute as well, at t 3> to- In 
Ref. [53[ it has been shown that these effects are negligible for 
t < 2to . In general the maximum admissible value of 4 at a given 
to depends on the underlying spectrum and on the basis of trial 
wavefunctions used. 
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and En- The correction factor arises from the finite di- 
mensionality and non-orthogonahty of the operator basis. 
IdeaUy one wih aim at a set of operators that dominantly 
couple to the first TV states and that are as orthogonal as 
possible to each other. The eigenvectors, up to the rota- 
tion and the change in the normalization of Eq. ([5]) , ap- 
proach their physical counterparts v" of Eq. Q too, with 
similar exponential corrections in Euclidean time 53]. 

From the eigenvalues we can also define effective energy 
levels, or, for p = 0, masses. 



f{t + a/2) = a"Mn 



A"(t,^o) 
A" -I- a, to) ' 



(8) 



that, for sufficiently large to and t > Iq, should exhibit 
plateaus which we then fit to a constant to obtain the 
masses to„ . 



S"1.4 



point 

naiTow 

wide 




2 4 6 8 10 

t/a 

FIG. 1. Effective masses of correlation functions between a 
point source and point, narrow and wide smeared sinks. 



C. Fermion field smearing 

The variational method is the central tool of our anal- 
ysis. It needs to be supplied with suitable building blocks 
in terms of operators, from which good approximations 
of the physical eigenstates can be obtained. The wave- 
functions of physical eigenstates will not be ultra-local 
objects and spatially extended interpolators need to be 
considered. We generate such extended operators by ap- 
plying Wuppertal smearing [s^ to a Fermion field tp: 



in) 



1 



1 + 6(5 



(n-l) 



±3 



(9) 



in the j-direction, x + aj, for instance a gauge link U^j- 
In our implementation we used APE smeared [55l | links 
for Ux,j, see Sec. Ill Dl below. Note that the Wuppertal 
smearing operator is gauge covariant. It transforms as a 
singlet under Oh, parity and charge transformations, it 
is Hermitian, translationally invariant and spin-diagonal. 



A 

or 

V 



- 0.01- 




0.001 



FIG. 2. Deviations of the average spatial plaquette from 
unity, against the number of APE smearing iterations Eq. (|16p 
on a lattice of ensemble for different a values. 

We can rewrite the above equation by defining a co- 
variant spatial lattice Laplacian, 



±3 



to obtain. 



1-1-6(5 



We introduce a fictitious time t — n„upAi, 

dm i'it + At) - m ^ k—w^m 



dt 



At 



At 



(10) 



(11) 



(12) 



where 



1 + 6S' 

The diffusion equation Eq. (jl2D is formally solved by. 



(13) 



m 



V'(O). 



(14) 



Starting from a (5-source "00; (0) = 5xO on a free configu- 
ration Uxj — 1 this results in a Gauss packet with the 
root mean square (rms) radius of V'^V'i 



n = 1 , . . . , TT-wup counts the iteration number and (5 > is 
a free parameter. The (arbitrary) normalization conven- 
tion is chosen to avoid numerical overflows for large iter- 
ation counts n^up- Ux,j is a gauge covariant transporter, 
connecting the lattice point x with its spatial neighbour 



^ = 3^kt/At = 3^ V^i;^ . (15) 

The diffusion speed is maximal for S ^ oo (k ^ 1/6) 
while the resulting wavefunction is more continuum-like 
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FIG. 3. The Wuppertal smearing function (riwup ~ 100, S — 0.3) with the original gauge hnks as paraUel transporters (left) 
and with APE smeared transporters (right). 



for (5 -> (fc — > 0). As a compromise we choose 6 = 0.3 
(fc « 1/9.3). 

By adjusting riwup we can control the overlap of our 
trial wavef unctions with the physical states. Using a 
point operator will lead to an effective mass with a signifi- 
cant curvature at small Euclidean times. A few iterations 
of smearing can help to flatten this out, suppressing the 
overlap with high excitations that have many nodes in 
their wavefunctions. Our strategy is to use an operator 
basis with a point operator that couples well to excited 
states, a narrow operator that couples well to the ground 
state and one operator that is somewhat wider. 

In Fig. [T] we display effective masses for the pseu- 
doscalar charmonium state, with a C75C point source and 
with a point as well as with smeared sinks, on ensem- 
ble (2), see Table m Note that since creation and destruc- 
tion operators differ, in the smeared cases the effective 
masses do not need to decrease monotonically. We em- 
ployed n„up — 20 and 80 smearing iterations for the nar- 
row and wide sinks, respectively. Note that we smeared 
both quark and antiquark fields so that the effective ra- 
dius of the charmonium creation operator is by a factor 
\/2 bigger than the expectation in Eq. ([T5|) . 



D. Gauge field smearing 

It was already suggested in Ref. 54] to replace the 
gauge links within Eq. ^ by other covariant transporters 
Ux,j, that depend on spatial links within the given times- 
lice. The ground state wavefunction is smooth and so we 
may wish to reduce the gauge field fluctuations as well, 
to enhance the overlap with low lying states. Follow- 
ing Ref. [13 we employ spatial APE smearing (551] to 
the gauge links that enter the Wuppertal smearing, it- 



eratively replacing a link by a linear combination of the 
link and the sum of the four surrounding spatial staples, 

-(/(«) _ r7("-l) I ^ \^ 7-r(»-l)7-r("-l) r7("-l)t 
'^x,i x,i "T " '-^x,j ^ x+aj,i^ x+ai,j ' 

C/i"^=^su(3)T4^.?. (16) 

a > is a weight factor and Pgu(3) A projects A onto U S 
SU(3) so that JieTr {U A^) is maximal. This procedure 
somewhat deviates from the definition of Ref. [l^l but 
also preserves gauge covariance. 

The spatial plaquette {Pa) measures the curvature of 
the gauge fields. Maximizing this means a smoother 
gauge background. In Fig. [51 1 — {Ps) is plotted against 
the number of APE smearing iterations on lattices of 
ensemble (see Table [I]) for three values of a. The 
approach to unity depends very little on the gauge con- 
figuration or on the gauge ensemble that we use. We 
decide to terminate the APE smearing after nape = 15 
iterations, using a = 2.5, as a compromise between gauge 
field smoothness and the computer time spent. 

APE smearing brings the links close to unity while pre- 
serving the gauge covariance of the Wuppertal smearing. 
This means that Eq. ([15]), which is valid for Ar ^ a 
on trivial gauge fields, is satisfied with good accuracy. 
In Fig. [3] we display a colour component after applying 
JT-wup = 100 smearing iterations to a (5-source in Coulomb 
gauge, without and with APE smearing. Indeed, the trial 
wavefunction looks smoother and moreover, we obtain 
the rms radius expected from Eq. (ITSt . 

Note that the APE smeared fields are only used to 
improve the operators but not to propagate the quarks. 
The inversions of the lattice Dirac operator are performed 
on the original gauge configurations. 



E. All-to-all propagators 

We first introduce the stochastic method to estimate 
all-to-ah propagators. We then describe the improve- 
ment methods that we employ, namely staggered spin 
partitioning (SSP) [4^ . the hopping parameter expansion 
(HPE) [56j and recursive noise subtraction (RNS) [il ]. 
We finally investigate the efficiency of combinations of 
these methods for a realistic example. Note that on 
top of this we also employ the truncated solver method 
(TSM) ^57, 58], see also Ref. that turns out to be 
beneficial even for masses as heavy as that of the charm 
quark. We restrict its use to light quark propagators 
though. 

1. Definitions and basics 
We denote the improved lattice Wilson-Dirac operator 

by, 

M =l-{l-nD). (17) 

This will depend on the quark mass through k. For each 
of the 12 (5-sources \0,a,a) at spacetime position 0, spin 
a and colour a we can compute solutions |s'^'"'°) of the 
linear systems, 

M|s°'"^°) = |0,a,a). (18) 
This defines the point-to-all propagator, 

5(^|0)X=5°'"'"(^,/5,&). (19) 

Due to translational invariance of expectation values, 
point-to-all propagators are often sufficient to calculate 
hadronic two-point Green functions. However, if one had 
all-to-all propagators at one's disposal, one would gain 
statistics from self-averaging over different source points. 
Moreover, some Wick contractions inevitably lead to dia- 
grams containing disconnected quark loops whose evalu- 
ations require more than a few source points. Solving the 
12 equations Eq. (fTO)) for all V lattice points (in our case, 
V = 131072 and 663552) instead of for a single = is 
prohibitive in terms of computer time and memory. 

However, we encounter statistical errors anyway from 
the path integral importance sampling in the calculation 
of expectation values. Hence it is sufficient to aim at an 
unbiased estimate, which can be obtained using stochas- 
tic methods [60, jil] . We introduce the following nota- 
tion, 

^ = ^''^=]^E^^ (20) 

and define random noise vectors \r]^), j — 1, . . . , N with 
components, 

r]'{x,a,a) = {x,a,a\if) e ^ (Z2 «) iZ2) . (21) 
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FIG. 4. Two dimensional sketch of a global noise source. For 
the propagator from x to y only the green line contributes to 
the signal, the black ones contribute to the noise. 



This complex Z2 noise has the properties, 

= O il/^/N) , (22) 

mif ^t + 0(l/VN) . (23) 

By solving, 

MW) = W) , (24) 

for |s^), j — 1,...,N, one can construct an unbiased 
estimate, see Eq. ([53]), 

M^' ■.= ]sM ^ M-\ (25) 
1 = - M-^ (1 - \T]){T]\j . (26) 

1 — \ri){ri\ is an off-diagonal matrix with entries of 
0{1/Vn). Hence the difference between the approxi- 
mation of Eq. (|25]) above and the exact result reduces 
like 1/Vn. When averaging over nconf gauge configu- 
ration the additional stochastic errors of an estimated 
observable reduce like l/\/Nnconf since the estimate is 
unbiased. One would like to achieve some sort of balance 
where this stochastic error becomes smaller than the un- 
avoidable gauge error cx 1 / y'riconf from averaging over 
the configurations. Since the ratio of these sources of er- 
rors is independent of riconf , once this is large enough for 
the central limit theorem to hold, this optimization can 
be performed on a small number of configurations. De- 
pending on the observable, a large number of estimates 
N may be required, unless the difference Eq. ((26)) can be 
reduced. Indeed, many methods of improving estimates 
exist, see, e.g. [5^ and references therein. 

Below we introduce the improvement methods that we 
use in this article. 
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FIG. 5. Two dimensional schematic sketch of standard spin FIG. 6. Two dimensional schematic sketch of odSSP (left) and 
partitioning. The numbers indicate the spinor component obdSSP (right). The numbers indicate the spinor component 
filled at the specific lattice site for set f (out of 4). filled at the specific lattice site for set 1 (out of 4). 



2. Staggered spin partitioning 

One source of large uncertainties of the naive estimate 
is that the noise source vectors have entries on aU lat- 
tice sites. The site, where the propagator ends, is sur- 
rounded by components of the source vector that may 
not contribute to the signal but that will contribute to 
the noise. To see this more clearly, consider the estima- 
tion of a propagator Eqs. (jlSp - (IT9t connecting the point 
X with the point y [see also Eq. (PB)) ]. 



SE{y\x)% = S{y\x)f^ 



(27) 



E 



1- |77)(77| (z\x 



/7Q : 



where only entries with either z^x^^^aoic^a give 
non- vanishing contributions to the noise sum, see Fig. U) 
Since signals decrease exponentially with Euclidean 
distances, 



\\S{y\z)\\ 



-\v-A/i 



(28) 



the source components located in the nearest neighbour- 
hood of y contribute most to the nois^ and thus it is 
desirable to reduce or to remove these terms entirely. 
Likewise, contributions that are off-diagonal in spin or 
colour at y should be avoided. A brute force way of 
achieving this is by "partitioning" [63 - l63 |. For the spe- 
cial case of spin partitio ning (SP) this is also known as 
the spin explicit method [63j. 

Partitioning amounts to decomposing TZ = V ® 



colour (g) spin into m subspaces TZf. TZ 



,iTZj. One 



can set the source vectors \riy) to zero, outside of the sub- 
space TZj and label the corresponding solutions as Is*^). 
The estimate of the all-to-all propagator is then given by 
the sum. 



(29) 



* This heuristic argument is over-simphstic since S is not a gauge 
invariant quantity. However, similar calculations can be per- 
formed for errors of physical observables, with the same result. 
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FIG. 7. The coupling strengths between the spinor compo- 
nents of standard (left), off-diagonal (centre) and off-block- 
diagonal (right) spin partitionings. Red indicates a strong 
(i.e. nearest neighbour) coupling, green a weak (i.e. next to 
nearest neighbour) coupling. 



Clearly, this results in an m-fold increase of the total 
number of solver applications. If the stochastic noise re- 
duction exceeds a factor l/y^ then this computational 
overhead is justified. 

Here we utilize spin and colour partitioning. So far 
within each spin partitioning set the same spinor compo- 
nent was dialled on every lattice site, see Fig.O Depend- 
ing on the observable, however, it may be favourable to 
alter the component to be filled within a specific set as 
a function of the spacetime position. For heavy quarks, 
the coupling between the upper and the lower two com- 
ponents of the Dirac spinor is small. One may exploit 
this by separating in spacetime components that couple 
strongly, only allowing for weakly coupled components to 
share a link. We call such non-trivial spin partitioning 
schemes staggered spin partitioning (SSP) [4J|. In Fig.[S] 
we sketch two SSP versions, off-diagonal SSP (odSSP) 
and off-block-diagonal SSP (obdSSP). Fig. [7] iUustrates 
the coupling strengths between the four spinor compo- 
nents. Red lines indicate nearest neighbours (strong cou- 
pling) and green lines next to nearest neighbours (weak 
coupling) . 

The obdSSP scheme should be particularly well suited 
to heavy quarks. However, this also depends on the 
discretization of the Dirac matrix and on the F- and 
derivative-structure of the creation operators. Hence pre- 
dicting the efficiency of a specific scheme is difficult. An 
object frequently appearing in this work is the pseu- 
doscalar loop Tr(75M~^). In our spinor representation 
75 is diagonal so that the naive picture presented above 
may apply. For other, non-diagonal F-structures differ- 
ent SSP schemes may be more effective. The picture be- 
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FIG. 8. Two dimensional sketch of the effect of HPE. y indi- 
cates the sink and x an arbitrary source site. One application 
of kD reduces the blue, two applications the green and three 
applications the yellow contributions to the noise. 



comes further obscured since within all the partitioning 
schemes there will be residual couplings between differ- 
ent colour components on the same site too. Fortunately, 
these terms are suppressed by the fact that S{y\y) will 
be quite colour-diagonal, in particular in the heavy quark 
limit. We also investigate combinations of (S)SP and 
colour partitioning. 



3. Hopping parameter expansion 

We have seen above that stochastic noise components 
that are close to the diagonal of the inverse Fermionic 
matrix are accompanied by larger amplitudes than 
terms that are far off the diagonal. Therefore the can- 
cellation of near-diagonal noise requires a comparatively 
larger number of estimates. The HPE noise subtrac- 
tion [56^ is based on the observation that one can expand, 
see Eq. P71) . 

oo k — 1 

= 2k ^(kD)* = 2k ^(kD)* + {nDfM-'^ , (30) 



i=0 



i=0 



where fc > 1. For distances between source and sink 
that are composed of more than k links, the first term 
on the right hand side does not contribute since D 
only connects nearest spacetime neighbours. Therefore, 
M~y = [{hD)'^ M~^]xy for sufficiently large source and 
sink separations. However their estimates will differ, 
M]^\y 7^ [{'^DY^E^]xy The variance of the latter esti- 
mate of M~y will be reduced since less noise terms con- 
tribute and in particular the dominant sources of noise 
have been removed. This was for instance exploited in 
Refs. [11,163. 

We illustrate the HPE technique in Fig. |8j one ap- 
plication of kD reduces the blue contributions, two ap- 



i - 
f 

i -0.01 - 




-0.02 -0.01 0.01 0.02 

Tr[ls><ri|Yj] 

FIG. 9. Scatter plot for the RNS, see Eq. ((ST)) for the pseu- 
doscalar loop. 



plications the green ones, three applications the yellow 
ones and so forth. However, note that unlike in the case 
of partitioning, these are not completely removed since 
they can propagate along a longer path to reach the sink, 
weakening their effect. It is self-evident that the HPE will 
be particularly effective for heavy quarks. 

Here we also study closed loops, i.e. x = y. Obviously, 
only even powers of D contribute to Tr(M~^r), where 
r e {1, 7m, CT^^^, 7^75,75}- We can write, Tr[M-^T) = 
K'^Tr (D'^M-ip) for k < fc^ax- For F = 75 and F = 7, 
for the clover action, /cmax = 2. Moreover, the lowest 
non-vanishing terms have been calculated analytically 
and can be computed and corrected for exactly (unbiased 
noise subtraction) [H, [66l - l68j . Here we do not attempt 
to do this but we restrict ourselves to fc = 2 instead. 

The HPE comes with very little computational over- 
head and, unlike in the case of partitioning, no additional 
solves are required. 



4- Recursive noise subtraction 

Within the RNS method 01 the off-diagonal terms of 
Eq. are estimated and subsequently corrected for, 

M-'=Wji\ + M-'{i-Wji\) 

^Wn\ + \sM{t-W)(ji\) , (31) 

in the hope to arrive at an improved estimate. The sec- 
ond term of the last line of the above equation involves 
additional inner products {rii\rij). For i 7^ j these fluctu- 
ate randomly but we sum over 12V 3> -/V^ such terms, 
where N is the number of estimates. Clearly, the pro- 
cedure can only work if this inner product is taken over 
a small er sub space. Therefore, we compute the random 
matrix |?7)(?7| only in the spin-colour subspace, setting all 
elements connecting different sites to zero. 
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TABLE IV. Gain factors of the noise reduction methods 
tested for the pseudoscalar disconnected correlator Eq. p2p 
at the charm quark mass, k denotes the number of kD apph- 
cations. 



k 





2 


no partitioning 


1 


2.89 


spin 


1.43 


6.32 


colour 


1.80 


5.06 


spin + colour 


2.52 


10.24 


odSSP 


2.30 


5.42 


obdSSP 


1.97 


7.16 


obdSSP + colour 


3.63 


11.80 


RNS 


1.87 


5.44 



In Fig. ini we display the correlation between the two 
terms of Eq. (PT|) for the pseudoscalar loop Tr(M~^75). 
Since the two quantities are anti-correlated, adding them 
together reduces the statistical error. So far we have 
assumed the coefficient of the second term of Eq. (|31l) to 
be unity. However, realizing that this term is an unbiased 
estimate of zero, one can generalize this method, e.g. by 
allowing for an adjustable coeffi cient. Moreover, difi^erent 
terms involving powers of 1 — |77)(77| and/or a different 
subspace where this matrix assumes non-trivial values 
could be implemented. However, we have not explored 
these possibilities any further. 



5. Comparison of methods 

We will apply the methods presented so far, to calcu- 
late disconnected quark loop contributions to the char- 
monium spectroscopy. An important example is the cor- 
relation of two zero-momentum projected disconnected 
loops, e.g., 

^ (Tr 1$),,75] Tr [{^M-'<^U^,] ) , (32) 

a:,y 

for the rjc mass where 1/4 ^ X4 + t. ^ = denotes 
a Wuppertal smearing function and the traces are over 
spin and colour. Note that the above traces are real, 
due to AI^ = 75Af75. We also remark that we used two 
separate sets of noise vectors to estimate the two traces, 
as one has to do. To compare the methods we choose 
t — a and n^up = 10 on ensemble 0. 

The improvements in terms of the real computer time 
spent to achieve the same stochastic error are displayed 
in Table IIVI for different combinations of methods. All 
overheads are included, except for the (negligible) cost 
of the two kD applications, k denotes the power of kD 
applied to the solution vector. The biggest net gain factor 
amounts to almost 12. Based on these numbers we decide 



TABLE V. Quark bilinears that we use (s^fe — \eijk\, also see 
Eq. (O)- 



name 


Oh repr. 




state 


operator 


TV 


Ai 


0-+ 


Vc 


75 


P 


Ti 


1— 


J/* 


7i 


bi 


Ti 


1+- 


he 


7i7i 


ao 


Ai 


0++ 


XcO 


1 


ai 


Ti 


1++ 


Xcl 


75 7i 


(P X V)t2 


T2 


2++ 


Xc2 


Sijfc7j Vfc 


(tt X D)t2 


T2 


2-+ 




7475 A 


(ai X V)t2 


T2 


2— 




75Sijfc7iVfc 


(P X D)a2 


A2 


3— 




7iA 


(61 X D)a2 


A2 


3+- 




74757^-0^ 


(ai X D)a2 


A2 


3++ 




757iA 


(ai X B)t2 


T2 


2+- 


exotic 


"/sSijk^jBk 


(61 X V)ti 


Ti 


1-+ 


exotic 


7475eijfc7iVfe 



to use the obdSSP and colour partitioning, together with 
the HPE for this type of diagram. 

III. THE SPECTRUM 

In this section we calculate the spectrum created by cc 
quark bilinears, neglecting charm quark annihilation and 
light quark creation diagrams. We first discuss our oper- 
ator basis and then the spectroscopy results. The varia- 
tional method also reveals information about the spatial 
structure of the underlying states. We will discuss this 
as well as the mixing between S- and D-wave operators 
in the J^'-^ = 1~~ vector channel. 



A. Extraction of masses 

The operators that we employ to create the charmo- 
nium states are based on Ref. [63. We restrict ourselves 
to the subset of these operators for which we are able 
to obtain meaningful signals. These are displayed in Ta- 
ble |Vl together with their irreducible lattice representa- 
tions and the corresponding lowest continuum spin as- 
signments, see Table HIl These assignments are of course 
not unique. For instance a radial Ti excitation can, in the 
continuum limit, very well correspond to J = 3 or J = 4, 
see Tabic Hm We label the operators by the names of the 
corresponding isovector mesons (which in nature are no 
charge eigenstates) since we are most familiar with these. 

Note that in the non-relativistic quark model, 

P={-)^+\ C=(-)^+^, (33) 

where S e {0, 1} and J G {L + S,L,\L - S\}. The 
states that cannot be accommodated in this way, namely 
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, 0"' , 1 ^, 2"' , 3 '", . . . are commonly referred to as 
"spin-exotic" states. As one can see from the 1 ^ ex- 
ample of Table |V] these exotic states do not need to be 
tetraquarks/molecules or hybrid mesons. With relativis- 
tic quarks, local bilinears are not restricted to ^ and 

1 anymore but e.g. 1^ can be created with L = 0. In 
this case the "exotic" 1 ^ state is merely a 1^ quark 
bilinear in a P-wave. We also remark that in QCD with 
finite quark masses L is not a good quantum number. 
However, it may still provide us with some guidance if the 
quarks are heavy. We will address this issue in Sec. IIII Dl 
below. 

In Table |Vl V represents a covariant spatial deriva- 
tive and D and B the symmetrized and anti-symmetrized 
combinations, 

Di = s^jkVjVk , Bi = f-ijkV-jVk , (34) 

with Sijk — l^ijkl and we sum over j and k. All operators 
containing a covariant derivative implicitly also include 
gluonic contributions but then any P- or D-wave will 
include derivatives and one would hardly call these hy- 
brid mesons. However, the _Bi-operators not only contain 
the vector potential but they are proportional to compo- 
nents of the field strength tensor itself. This is as close 
to a valence gluon as one can get. A strong coupling of a 
physical state to this operator may then indicate a large 
hybrid meson content. The 1 ^ charmonium is consid- 
ered a prime hybrid candidate. However, we find the 
operators SijkJjBk (Ti representation, not listed in the 
Table) to be very noisy, with no compelling evidence of 
a coupling to the ground state created by 6i x V. 

For each operator listed in the Table we construct 
a three by three cross-correlator matrix, see Sec. IIIB[ 
with different smearing levels, see Sec. Ill CI We apply 
the same smearing to quark and antiquark. The smear- 
ing parameters have been optimized for several states 
as described in Sec. Ill CI so that point-smeared effective 
masses are relatively constant for the narrow operator 
and approach their asymptotic values from below for the 
wide operator. We only consider the two lowest lying 
masses reliable since the highest state contained in the 
basis may be polluted by even higher excitations. In 
Ref. (t^ a similar approach was used to calculate the 
spectra of excited states in the quenched approximation. 



B. Discussion of the results 

We display effective masses for ensemble (see Ta- 
ble m , obtained after diagonalizing the correlation ma- 
trices at the normalization time to, see Eq. (O, for some 
of the channels in Fig. [TOl Only the lowest two masses 
are fitted and the fit ranges are indicated by the blue 
lines. The extracted mass values, together with these fit 
ranges, are displayed in Table IVTl In this table we also 
assign the lowest possible continuum spin to each chan- 
nel. Note however, that radial excitations of both J = 1 



TABLE VI. Fitted masses obtained on ensemble for the 
first two eigenvalues in each channel. The normahzation time 
to and the corresponding fit ranges are also given. The errors 
are only statistical and we give the lowest continuum J^^ 
from which the lattice representation can be subduced. 



operator 


,7^"^ 


to/a 


mi/McV 


range 


m2/MeV 


range 


TT 


0" + 


1 


2993 (4) 


5-12 


3645 (19) 


1-8 


P 


1"~ 


1 


3070 (6) 


7-12 


3699 (24) 


1-7 


bi 


1 + - 


2 


3457 (22) 


2-7 


4060 (65) 


1-5 


ao 


0++ 


2 


3381 (19) 


4-12 


3996 (48) 


1-5 


ai 


1++ 


2 


3462 (20) 


3-11 


4011 (52) 


1-5 


(P X V)t2 


2++ 


1 


3471 (19) 


1-6 


3917 (46) 


1-6 


(tt X D)t2 


2- + 


1 


3756 (32) 


1-9 


3995(141) 


1-6 


(ai X V)t2 


2"" 


2 


3706 (27) 


1-10 


4076 (83) 


1-6 


(P X D)a2 


3"" 


1 


3782 (35) 


1-8 


4815 (92) 


1-6 


(61 X D)a2 


3+- 


1 


3995 (50) 


2-6 


5365 (76) 


1-3 


(ai X D)a2 


3++ 


2 


3993 (54) 


1-5 


5008(287) 


1-4 


(bi X V)ti 


1- + 


1 


4154 (54) 


1-5 


4297(181) 


1-4 


(al X B)t2 


2+- 


1 


4614(220) 


1-9 


4643(254) 


1-8 



{— Ti) and J ~ 2 {— T2) states could in principle corre- 
spond to J = 3. In particular, this possibility cannot be 
excluded for the excitations of the Ti states that we have 
labelled as 1^ and 1++ and for the T2 state labelled as 
2++. 

The ground states and first excitations of the stan- 
dard S- and P-wave states rjc, J/^*, he and XcJ exhibit 
good signals and stable plateaus. The effective masses of 
higher spin states are naturally noisier and thus compli- 
cate the fits. 

A particularly interesting channel is the 1 '". Al- 
though this is a spin-exotic state it can be well accessed 
by the 5i x V operator that does not contain an explicit 
chromo-magnetic field. However, V contains a link vari- 
able and may allow for a gluonic excitation. The quality 
of the effective masses arising from the &i x V operator 
is not high but the fit yields reasonable errors. The two 
lowest lying states are very close. In fact, within their 
errors the effective masses are overlapping so that in the 
statistical analysis it was necessary to sort the jackknifes 
according to the proximity of the eigenvectors to the ones 
obtained on the original ensemble. This may hint at a 
hybrid nature of this channel. Static hybrid potentials 
are repulsive at short distances so that, within a poten- 
tial model, we ma y expe ct smaller energy gaps between 
radial excitations |7ll - l73| . 

The computed spectrum is displayed in Fig. 1111 to- 
gether with the experimental valuefO We have used the 
spin-averaged IS* mass Eq. ^ to fix the charm quark 
mass. The other states are predictions. Note, however, 
that we underestimated rn^^ by about 15 MeV. The main 
effect of this is that all the predictions should be shifted 



^ In some cases their J^^ assignment is still under debate. For 
instance for the X(3872) that we list as a 1++ state a 2 ^ as- 
signment is not ruled out. 
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FIG. 10. Effective masses for some operators of Table IVl obtained on ensemble 0. The fit ranges and errors are indicated by 
horizontal lines. The to values (in lattice units) refer to the respective normalization times Eq. 



up by 15 MeV. Keeping this in mind we observe all spin- 
averaged states below threshold coming out in qualitative 
agreement with the experimental data [HO]- We indicate 
the experimental DD and DD open charm thresholds 
as horizontal lines. Negative parity states cannot decay 
into DD and the DD threshold is believed to play a 
more prominent role in the decay of hybrid mesons than 
DD*. 

The spin-averaged IP-lS* splitting is underestimated, 
relative to experiment while the 2S-1S splitting comes 
out right. There may be issues with the scale setting, 
due to the unrealistic sea quark content. We have also 
remarked in Sec. Ill Al above that there are reasons to 
believe ^8J that the lattice spacing should be set to « 
1.81 GeV rather than to the a'^ sa lJ2_GeV that we 
used. This change would bring the IP- IS* splitting in 
line with experiment but result in an overestimated 25* 
mass. This in turn could then be due to a combination of 
finite size effects and interferences with the DD threshold 
in nature that do not occur on our ensemble, due to the 
heavy sea quark mass. 

From pNRQCD one would, to leading non-trivial or- 



der in l/rric, expect the S'-wave finestructure to be de- 
termined by the matrix element [t^, [zE| j 



where to leading order in perturbative QCD, 



(35) 



(36) 



^ff is the non-relativistic charmonium wavefunction, as 
the strong coupling parameter and cf = 1 + C(as) is a 
matching coefficient that has only recently been deter- 
mined in lattice schemes ;76]. 

This illustrates the very short-distance nature of the 
S'-wave finestructure that should be affected significantly 
both by lattice spacing effects and by differences in the 
running of the coupling depending on the sea quark con- 
tent. It is also very sensitive to the charm quark mass. 
Reducing this by 5 % would increase the splitting by 
10 %. So it is not surprising that we underestimate the 
experimental number of 117 MeV for the IS* finestructure 
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FIG. 11. The predicted spectrum, together with the experimental values on ensemble 0, see Table|Tl 



splittmg. We obtam Amis = mj/^, — m,,^ = 77(2) MeV 
on ensemble 0, Ami5 = 88(4) MeV on ensemble (2) 
and Atoi5 = 130(9) MeV on ensemble (J). Indeed, with 
lighter sea quark masses this value increases. 

For the 25 finestructure splitting we obtain Am25 = 
™*(2S) —iT^v'a ~ 54(6) MeV on ensemble and Am2s = 
56(8) MeV on ensembleH (2), in agreement with the ex- 
perimental value of 49(4) MeV. In view of the disagree- 
ment of the IS splitting this is quite surprising since one 
would have expected a lot of the systematics to cancel 
from the ratio of the 25" hyperfine splitting over the IS* 
splitting, see Eq. ([35]) . We may therefore wonder whether 
either the physical rjc or the ^'(25) states are unusually 
low, due to contributions from quark line disconnected 
diagrams. In the first case our neglection of cc annihila- 
tion diagrams may be relevant while in the second case 



On ensemble (J) where the radial excitations are seriously af- 
fected by the finite volume we get Am2s = 177(66) MeV. 



omitting qq creation (and the use of unphysically heavy 
light quark masses) would be the dominant efFect(s), see 
Sees. IIVI and |V] below, respectively. 

We remark that we also underestimate the P-wave 
finestructure. This is expected too and again mostly due 
to lattice spacing effects and an unrealistic sea quark con- 
tent. We also notice that in our approximation where the 
open charm thresholds are much higher than in nature 
the Z(3934) (recently renamed into Xc2(2P) [13]) niay in- 
deed be associated with the x'c2 state while the X(3872) 
certainly is lighter than one would have expected from 
an excited P-wave. However, in the first case we cannot 
exclude the possibility that we have misidentified a 3+^ 
state as 2"'""'", in particular since this comes out lighter 
than the other two x'c niultiplet masses. Finally, the 
proximity of the two 1 ^ states as well as of the two 
2"' states may indicate a hybrid nature of these spin- 
exotic charmonia. We have not detected such indications 
in any of the other channels. With the exception of the 
A2 (J = 3), in these cases the radial excitations are lower 
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in energy than these spin-exotic states. 



"Wavefunctions" 



In Sec. Ill B I we have introduced couplings between an 
operator Oi, i = 1, . . . , N , and a physical state vf = 
{0\Oi\n). These will be approximated, up to a rotation 
and normalization, see below, by the ip'^{t, to) of Eq. ([5]). 
We employ the normalization, IV'"(^j^o)P = 1- In the 
pseudoscalar channel our operators read. 



Ot = Yl c(y)*»(y ~ x)75c(x) 



(37) 



Here denotes the square of the Wuppertal smearing 
operator, Eq. ©, since we apply this to quark and anti- 
quark fields. $i is translationally invariant and will only 
depend on the difference y — x. 

We employ a iV = 4 dimensional trial basis consist- 
ing of TT-wup = 0, 5, 10 and 40 iterations. This means that 
the contain twice these numbers of iterations. Folding 
these smearing functions with the asymptotic eigenvec- 
tors results in a new smearing function. 



where 



(38) 



(39) 



U G SO{N) diagonalizes C{to) and dj > are the square 
roots of its eigenvalues. 



C-i/2(io)C(t)C-i/2(to) 



(40) 



see Eq. ([S]). In particular this means that the operators 
constructed from <I>j, Oj [see Eq. ([57| ] create states that 
are orthonormal ai[3 t — to'- {Oi{to)Oj{0)) — Sij. Also 

note that if we neglect the coupling of the operator Oi to 
states with energies bigger than En, di cx exp(— £"^^0/2). 

If we had perfect overlap with the respective physi- 
cal states then we could choose = t = Q. In this 
case, in the non-relativistic limit, where we do not en- 
counter particle-antiparticle creation, one may identify 
|vl/"t(x)vl/"(x)| with the quantum mechanical probabil- 
ity density. On a qualitative level, we may still think of 
4'"(x) as the wavefunction of the nth state. The used 
ensemble is unfortunately too coarse to resolve the 



If corrections from truncating the basis can be neglected for 
to = (which is unlikely) then, at this ip, t/"" — > v" at large 
times t. 



node structure of the gauge invariant |\I/"t\[/n|^ However, 
one can also plot a diagonal colour component of (x) , 
after fixing to Coulomb gauge. In fact the APE smeared 
gauge links are so close to unity that it is hard to re- 
solve the differences between Coulomb gauge fixing and 
a non-gauge fixed component from a plot. 

In Fig. [TH we show a two dimensional cross section of 
one colour component of the normalized wavefunctions 
^'"(x), n = 1, 2, 3. In spite of the small basis and lattice 
volume the node structure is consistent with the 15', 25* 
and 35* assignments, with no visible pollution from higher 
Fock states or Z?-waves. For the 15 "wavefunction" we 
obtain an rms radius, Ar = (r)rms = t'^I^'P ~ 

0.39 fm. This compares reasonably well with the infinite 
volume continuum potential model expectation of about 
0.4 fm 14911. 



D. Mixing in the vector channel 

Due to its direct production in electron-positron an- 
nihilation, the vector channel is rich in experimentally 
confirmed resonances, as can be seen in Fig. [TT] Of great 
interest is the inner structure of these states, in particular 
of the 4* (25) and the 5* (3770) states, which have a mass 
difference of only about 90 MeV and both are close to the 
DD open charm threshold. While J/\l/ is dominated by 
15 quark-antiquark configurations, its excitations may 
exhibit a more complex structure. 

As the name suggests, the \l/(25) is thought to be a 
radial excitation. Since \E'(3770) is so close in mass, it is 
very improbable that it is excited in a further, higher ra- 
dial vibration mode. One possibility which we investigate 
here is an orbital excitation where the quark-antiquark 
pair is in a relative £)-wave. We start from an operator 
basis consisting of three 5-wave and two £)-wave inter- 
polators, 

(c7ic)o , (c7ic)2o , (c7ic)8o , (41) 
{cSijkljDkc)o , (cSijk^jDkCjso , 

where Dk is defined in Eq. (j34l) and Syfe = \^ijk\- The 
subscripts denote the numbers of smearing iterations, 
both for the quark and antiquark. Initially we planned 
to include hybrid operators like cj^BiC into our basis but 
unfortunately these provided very poor signals through- 
out, independent of the smearing levels. 

This mixing analysis is performed on ensemble (2), see 
Table |T1 with to = 3a. We display the lowest four ef- 
fective masses in Fig. [T31 Indeed, the second and third 
eigenvalues lie very close. The fourth eigenvalue may be 
identified with the ^"(4040). 

The eigenvector components reveal the overlaps be- 
tween the trial operators and the physical eigenstates. 
The correlation matrix that is real in our case has the 
normalization ambiguity Cij{t) 1— CiCjCijit), where i — 
1,. . . ,N and N is the dimension of the operator basis. 
One convenient choice is Cii(O) = 1. The orthogonal 
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FIG. 12. The IS, 25' and 3S pseudoscalar "wavefunctions" 
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FIG. 13. Effective masses of the four lowest lying states in 
the vector channel on ensemble &)■ 
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FIG. 14. Components of the first 1 
Eq. gH). 



eigenvector, see 



transformation U that diagonalizes the correlation ma- 
trix at the time to and the eigenvalues at this time of 
Eq. (PP]) depend on this initial normalization of C(i). Fol- 
lowing the discussion of Sees. Ill Bl and IIII Cl see Eqs. ([5]) 
and (j40p . we can define effective overlaps, 



-1-1/2 



3 ■' 



(42) 



that in the limit to ^ 0, t ^ will approach = 
(O|0i|7i), up to an overall factor. The effective over- 
laps do not depend on the normalization choices e.i of 
Cij{t). Our normalization X^i kr(^' ^o)P = 1 also im- 
plies J2n \^?i^'^o)\'^ — 1 which is equivalent to ignoring 
any effects of higher lying states. In Figs. [TH ITSlandfTBl 
we display the first three v^{t,to = 3a). 

As one would expect the ground state, the J/'^, does 
not receive any contributions from the two D-wave op- 
erators. Interestingly, the second eigenstate, that is en- 
ergetically very close to the third one (in fact for most 
t-values we can only differentiate between these states 
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FIG. 15. Components of the second 1 eigenvector. 
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FIG. 16. Components of the third 1 eigenvector. 



by tracing their eigenvector components), does not "see" 
any D-wave operators either. We remark that at the t- 
values where these second and third energy eigenvalues 
differ from each other we encounter more "mobihty" of 
the eigenvector components, see Figs.[T21[Tn]and[Tni Note 
the relative sign change in the case of the first excitation 
between the local/narrow and the wide operators, result- 
ing in a node of the spatial wavefunction, similar to the 
25 state of Fig. [HI This strongly suggests a ^"(25) as- 
signment for this state. Conversely, the third eigenvalue 
only couples to the wide smeared D-wave operator, which 
obviously leaves it as a candidate for the ^'(3770). These 
results compare reasonably well with the ones of Ref . [t^] ■ 
We conclude that the charm quark is sufficiently heavy 
for S- and D-waves to undergo only very mild mixing. 
So, at least for charmonia of masses below 3.8 GeV, it 
is meaningful to label states by their angular momenta. 
However, we have not yet considered the effect of open 
charm thresholds. We will address this question in Sec.lVl 
below. 




FIG. 17. The lowest order perturbative QCD graph responsi- 
ble for the mixing between C75C and q-ysq states. The red lines 
correspond to charm quarks, the black ones to light quarks, 
the twiddly ones to gluons. 



IV. MIXING BETWEEN THE rye AND THE 
LIGHT T7 MESON 

Charmonia are fiavour-singlet states, however, so far 
we have neglected the charm quark annihilation diagram 
that arises from Wick contracting the correlation func- 
tion, ([crc](t) [crc]^(0)). The inclusion of quark line dis- 
connected diagrams may affect charmonium masses. In 
particular the proximity of the mass of the rjc meson to 
that of the pseudoscalar glueball which may propagate as 
an intermediate state after cc annihilation may have some 
effect [77]. This glueball mass was consistently deter- 
mined in simulations of pure gauge theories on isotropic 
lattices to b e 1781 (2630 ± 290) MeV and on anisotropic 
lattices [7^^(2637 ± 26) MeV. We also know that 
the light quark analogue of the rjc, the 77' meson is much 
heavier than the light octet pseudoscalar mesons, due to 
the Ua(1) anomaly. Naturally, chiral symmetry will not 
play a prominent role for charm quarks. However, this 
does not exclude a remnant effect of the vacuum topology 
that may lift the rjc mass by a few MeV. 

First calculations of the disconnected contribution 
both with np — 2 sea quarks and in the quenched ap- 
proximation are consistent with no mass shift of the rjc 
mass [80, 81], within statistical errors of about 20 MeV. 
More recently, in the quenched approximation the r]c 
mass was estimated to increase by 1-2 MeV [s^] due to 
disconnected diagrams and, including sea quarks, this ef- 
fect may become 1-4 MeV [82| . 

Naturally, when sea quarks are included both physical 
eigenstates, the rjc and the 77', will contain light as well 
as charm valence quarks. The charm-anticharm compo- 
nent will be dominant within the r/c while the r/ will 
almost exclusively contain light quarks. In our case we 
employ np = 2 sea quarks and therefore an isosinglet 77 
state assumes the role of the 77'. In addition we have an 
isovector tt triplet, instead of the octet. When the dis- 
connected quark loop is included then, at large Euclidean 
times, the ijc state will decay to the ground state in the 
J^^ = ^ channel, which is the rj meson. The phys- 
ical rjc will only appear within the tower of excitations 
of this ground state. Following Ref. [i^l we call this ef- 
fect "implicit" mixing: the C75C state already intrinsically 



We converted the numbers into units of rg 
the uncertainty in this scale setting. 



: 0.467 fm and ignore 
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contains a 975(7 contribution. However, one would expect 
the coupling of the C75C creation operator to this state to 
be very weak. Otherwise charmonia would not be stable 
in nature either. This means that we can treat this as 
a perturbation. We decompose the physical Hamiltonian 
H = Hq + \Hi + • • • , into a part i/o witlj^ cc and qq 
eigcnstates, without pair creation. The small perturba- 
tion A-ffi is then responsible for the mixing. Neglecting 
radial, gluonic or multiquark excitations, the physical rjc 
wavefunction of this two-state system reads, to first order 
in the small parameter A, 



\'nc 



I- \ , X {m\Hi\cc) 

\cc) + A— — -^m) 



E- — R- 



(43) 



While we do not know the functional form of Hi or of 
the unperturbed wavefunctions, we can evaluate all the 
relevant matrix elements on the lattice. Fig. [T7] depicts 
the graph responsible for this mixing to lowest order in 
perturbative QCD. 



^^^^^^^^^^^^ - ^^"^"^ ^^^^ 


2 ^^^^^^ ^^^^^ 


2 '(^^^ 





FIG. 18. Correlation matrix for the mixing between states 
created by C75C and by qysq operators. Red lines represent 
charm quarks, black lines light quarks. 

Obviously, the mixing will depend on the light quark 
mass value m^. With decreasing the denominator 
of Eq. (|43)) will become larger. However, we would also 
expect the matrix element in the numerator to increase 
since the probability of creating a light quark-antiquark 
pair may increase with decreasing quark mass. Thus, 
ideally one would realize several light quark masses to 
clarify this issue. 

Similar to our investigation of mixing in the vector 
channel of Sec. HHP I we also calculate a correlation ma- 
trix here. We choose the basis states, 



(c75c)o , 

(9759)0 : 



(c75c)io 

(9759)5 



(C75C)80 • 

(9759)40 • 



(44) 



where the subscripts denote the number of Wuppertal 
smearing iterations. The diagonalization of this matrix 
at large times will not only allow us to extract the energy 
levels but it will also provide us with qualitative informa- 
tion on the charm and light quark content of the physical 




FIG. 19. The lowest order perturbative QCD graph respon- 
sible for implicit mixing between states created by C75C and 
by 5755 operators. 



states. In Fig. [18] we sketch the structure of the mixing 
matrix between the cc and qq sectors, omitting the dif- 
ferent smearing levels. Red lines represent charm quark 
propagators and blue lines light quark propagators. The 
prcfactors are due to the np — 2 mass degenerate sea 
quark flavours. The upper left corner contains the cc sec- 
tor, the lower right corner the qq sector. The off-diagonal 
elements quantify the mixing. 

Note that already in the pure charmonium sector "im- 
plicit" mixing will occur, due to intermediate light quark 
loops in the disconnected part, see Fig. I19[ or even inter- 
mediate glueball states. If the explicit mixing encoded 
in the 0(A) off-diagonal elements of the mixing matrix 
is small then we will not be able to resolve the O(A^) 
decay of cc states into states dominated by qq within any 
sensible Euclidean time distances. This further justifies 
and motivates our mixing matrix approach. 

Our strategy is as follows. We first determine the 
eigenvalues and eigenstates of the three by three sub- 
matrices within each of the flavour sectors, separately, 
in order to obtain an "unperturbed" approximation to 
the spectrum. We will then compare the spectrum and 
the eigenvector components of this reference point to the 
situation with the mixing elements switched on. 

The all-to-all propagator estimates for both, the charm 
and light disconnected loops have been improved by the 
HPE, obdSSP and colour partitioning, see Sec. Ill El For 
the hght quark propagators, in addition the TSM [57, 5^ 
with Ut = 25 has been applied. We analyse ensemble 0, 
see Table U with a pseudoscalar mass of about 1 GeV. 
At this heavy mass we find an 77- tt mass splitting of only 
52(13) MeV. Within the diagonal three by three subma- 
trices we find the disconnected charmonium loops to be 
very noisy. Since the statistical errors are bigger than 
the expected splitting of a few MeV we ignore these con- 
tributions. If we were to detect significant off-diagonal 
contributions to the full correlation matrix then of course 
we would have to revisit this issue at a later stage. 

The masses of the light 77 meson and of its first ra- 
dial excitatiorF^I 77' can be extracted from the largest two 
eigenvalues of the submatrix containing the light quark 
creation and annihilation operators while the 77c and 77^ 



^ We omit the F structure for convenience. 



This should not be confused with the pseudoscalar flavour-singlet 
meson in the rip =2 + 1 theory. 
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FIG. 20. Effective masses from the eigenvalues of the full 
matrix. As a reference point the effective masses from the 
(unmixed) submatrices are plotted too (black squares). 
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FIG. 22. Eigenvector components of 77c in the full basis. 



TABLE VII. Fitted eigenvector components of the rj and 77c 
states from the diagonalization of the full matrix. 
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FIG. 21. Eigenvector components of 77 in the full basis. 



masses can be approximated from the cc sector. We find 
a diagonalization of the full six by six matrix to be numer- 
ically unstable and hence restrict ourselves to the basis 
of the states (0750)10, (0750)80, (9759)5 and (9759)40 for 
our subsequent full-fledged mixing analysis. 

In Fig. [20] we display the effective masses of the ry, ry', 
77c and 77^ states, with the off-diagonal matrix elements 
switched off (squares). These are compared to the low- 
est three effective masses obtained from the four by four 
matrix with the mixing switched on, as explained above. 
We shift the latter effective masses slightly to the right. 
No relevant deviations can be seen. 

The mixing can be studied in more detail by investigat- 
ing the respective effective eigenvector components that 
are defined in Eq. where U diagonalizes C(to) that 
has the eigenvalues df. We display these effective over- 
laps for the ground state 77 meson in Fig. [21] and for the 
77c in Fig. [21] The fitted components are displayed in Ta- 



ble I VIII Indeed, the 77 does not receive any statistically 
relevant cc contribution and vice versa. The summed 
probability to find the 77 meson in a state that can be 





(cc) 10 


(CC)80 


{qq)s 


(95)40 


V 


-0.017(37) 


0.009(63) 


-0.806(1) 


0.591(9) 




0.333(30) 


0.943(11) 


0.000(41) 


0.021(47) 



created by the co operators amounts to (4 ±25) • 10~^ and 
to find the 77c in a state created by qq to (4 ± 22) • 10~^. 
These tiny upper limits on possible mixing effects also 
render a relevant coupling of the 77c state to the pseu- 
doscalar glueball extremely unlikely since this glueball 
can appear as an intermediate state in diagrams of the 
type depicted in Fig. [17] 

Obviously, the energy shift from explicitly admitting 
CO annihilation and light quark creation cannot signifi- 
cantly differ from zero either. We find a mass difference, 
^mixod _ ^unmixed ^ i^2A) McV. After this demonstra- 
tion of the feasibility of such studies we wish to further 
reduce the statistical errors and to vary the light quark 
mass in the near future. We address contributions from 
higher Fock states that may be relevant, e.g. for the 77^ 
state in the following section. 



V. MIXING BETWEEN cc AND DD 
MOLECULAR OR TETRAQUARK STATES 

Charmonia can decay into pairs of (excited) D and 
D mesons if their masses are above the allowed de- 
cay thresholds. Charmonia near these thresholds may 
however also contain significant Fock admixtures of DD 
molecules, see e.g. Ref. [83|, or of cqqc tetraquarks. We 
will study these effects in three different J^*^ channels, 
^, 1^^ and 1++. The first two channels are inter- 
esting with respect to the experimental overpopulation 
of the vector channel and the fact that the '^{2S)-r]'^ 



TABLE VIII. r structures of meson and molecule interpolat- 
ing fields, see Eqs. (|45l) and (|46)) . 
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finestructure splitting is very small, compared to that 
of the ground states, see the discussion of Sees. IIII Bl and 
IIVI above. The l"'"^ is phenomenologically relevant to 
disentangle the nature of the X(3872) state. 

We will address the question of higher Fock state con- 
tributions to the spectrum by creating and destroying 
states employing both, the cc operators corresponding 
to r^c, J/^ and Xci charmonia as well as the four-quark 
operators corresponding to DiD* in J ~ 0, DiD and 
D*D molecules, respectively. The analysis method is ex- 
actly the same as outlined in Sec. IIVI However, in the 
present situation the dependence of the mixing strength 
of Eq. (j43p (replacing cc i— > cqqc, qq i— > cc) on the light 
quark mass is evident: again, with decreasing light quark 
masses the numerator is likely to increase, however, the 
denominator will decrease as the energy gaps between 
open charm states and the first radial charmonium exci- 
tations become smaller. Therefore, we analyse ensemble 
(3) (see Table HI) that, with a light pseudoscalar mass of 
about 280 MeV, is closest to the physical point. Note that 
with the product mp^L sa 2.6 this L 1.84 fm lattice is 
quite small so that in particular for radial excitations we 
may expect significant finite size effects. 

We start from a six dimensional operator basis con- 
taining three cc and three molecular interpolators, differ- 
ing by their Wuppertal smearing levels. We label these 
as p(oint), n(arrow) and w(ide). The generic form of 
the meson operators centred at the spacetime position x 
reads, 

= (cFmc)^ , (45) 

where for readability we omit the smearing functions. 
Within the molecular operators we allow for a spatial 
separation r, 

YAv) [{qT\c)^{cT\q),+, 

+ {-Y {cV\q)x[qT\c)^+,] . (46) 

The r structures and s e Nq values for the states of 
interest are displayed in Table IVIIII see also Ref . [sj] . 
Note that in this exploratory study we restrict ourselves 
to operators with molecular contractions and ignore the 
possibility of arranging the quarks into tetraquark-like cq 
and qc diquark-antidiquark structures. 

In Fig. 1231 we sketch the structure of the mixing matrix. 
The different smearing levels are again omitted for clarity. 
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FIG. 23. Correlation matrix for the mixing of charmonia with 
DD molecules. Red lines represent charm quark propagators, 
blue lines light quarks. 



TABLE IX. Mass spectrum in MeV, neglecting the mixing 
between the two- and four-quark sectors. 





0-+ 


1"" 


1++ 


cc ground state 


3045 (8) 


3175(10) 


3660 (40) 


cc first excitation 


3991(46) 


4168(67) 


5043(120) 




4848(39) 


4715(39) 


4152 (38) 


-D(i)-D'' ' molecule 


4822(23) 


4712(23) 


4064 (28) 



Red lines represent charm quark propagators and blue 
lines light quark propagators. The prefactors are due to 
the n-p = 2 mass degenerate light sea quark flavours. The 
upper left corner contains the cc, the lower right corner 
the molecular sector. The off-diagonal elements encode 
explicit mixing. In our calculation we omit the charm 
annihilation diagrams of the second lines within each of 
the correlation matrix sectors; based on our experience 
of Sec. IIVI above we deem these numerically irrelevant. 
A similar matrix was constructed, e.g. in Refs. [4l|,[8^ in 
order to investigate the p meson decay width and light 
tetraquark states, respectively. Note that all the depicted 
diagrams that include light quark propagators include 
more than one quark line contraction since for r ^ 
Eq. pS)) contains two terms. 

The spatial separation within the molecular operators 
was tuned to maximize the correlation function of the 
molecular sector. This led us to employ the on- axis sep- 
aration r = 4a 0.3 fm. After some experiments we 
decided to employ point-to-all propagators for all dia- 
grams, with the exception of the top right diagram within 
the molecule-to-molecule sector, see Fig. [231 This neces- 
sitates to implement a light all-to-all propagator at the 
sink. For this purpose we generated the equivalent of 100 
complex Z2 stochastic estimates, employing the obdSSP, 
HPE and TSM methods. 

We first diagonalize the submatrices separately within 
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FIG. 24. Mass spectrum from separately diagonalizing the 
submatrices within the different sectors. The right points for 
the DD states correspond to the sums of non-interacting D 
mesons, the left points to interacting I — D and D mesons. 



the cc and molecular sectors to obtain a reference spec- 
trum. This provides us with up to four reliable eigen- 
values, two within each sector. However, the excited 
molecular channels are quite noisy so that in these cases 
we are only able to extract acceptable plateaus for the 
ground states. The remaining three states within each 
of the J^'~' channels are plotted in Fig. [M] Next to the 
isospirF^ / = molecular masses we also display the 
sums of the masses of the corresponding individual D and 
anti-Z? mesons. The resulting masses are also displayed 
in Table IIXI where the errors are statistical only. Due to 
the finite volume the radially excited S'-wave states are 
significantly higher than the corresponding experimental 
masses and the excited P-wave suffers even more from 
this effect, being by almost 1 GeV heavier than the cor- 
responding D* D molecule that is already heavier than in 
the real world due to the unphysically large light quark 
mass. 

Within errors of about 30 MeV we do not find any 
significant mass differences between molecules and their 
open charm constituent mesons in the pseudoscalar and 
vector channels. Of particular interest is the mass of the 
1++ molecule that is by almost 200 MeV heavier than 
the X(3872). However, this can easily be attributed 
to the light quark mass since the light pseudoscalar is 
still by 140 MeV heavier than the physical one. We 
find a significant binding of this axialvector molecule, 
m 



D'D ' 



{niD' +mD) = 88(26) MeV. There wih be some 
volume and light quark mass dependence of this value 
that needs to be studied. Note that this binding en- 
ergy is much bigger than the mass differences between 
the experimental X{3872) of a fraction of an MeV with 
respect to electrically neutral open charm states and of 



O J/V|/ 

D,D, 



8000 
7000 
6000 
5000 > 
4000 g 
3000 
2000 
1000 



5 6 7 8 9 
t/a 



10 11 12 13 14 



FIG. 25. Effective masses from the eigenvalues of the full 
matrix in the 1~~ channel. As a reference point the effective 
masses from the submatrices are plotted too (black symbols). 
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Note that the 1 = 1 channel has been studied by Liu [§6 



FIG. 26. Eigenvector components of J/^ in the full basis. 



a few MeV with respect to charged D and D* mesons. 
However, on a qualitative level the increased attraction 
deserves some attention. 

We base our fully fledged mixing study on a four by 
four submatrix with the operator basis Mp, Mn,Yp and 
Yn- The normalization time is to = 2a for all channels. 
We start the discussion with the vector state and display 
the resulting lowest two effective masses, together with 
the unmixed reference masses in Fig. 1251 The correspond- 
ing effective overlaps Eq. (02]) arc plotted in Fig. [25] for 
the J/^- and in Fig. [22] for the DiD molecule. The J/* 
receives the dominant contributions from the cc sector. 
However, the molecular configurations contribute signif- 
icantly too. In contrast, the DiD state only contains a 
small (but non-vanishing) cc admixture. This is very sim- 
ilar to the observation of Ref. [12] that the ground state 
potential between two static sources Q and Q receives 
a significant light quark contribution also for distances 
much smaller than the string breaking distance while its 
QqqQ excitation contains almost no QQ component. We 
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FIG. 27. Eigenvector components of DDi in the full basis. 



FIG. 29. Eigenvector components of xa in the full basis. 
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FIG. 28. Elfective masses of the eigenvalues of the full matrix 
in the 1^"*" channel. As a reference point the effective masses 
from the submatrices are plotted too (black symbols). 



ip 

0.8 - 

0.6 i 5 

0.4 - 
0.2 - 

- 

-0.2 -J ^ 
-0.4 - 
-0.6 - 
-0.8 ^ 

-i4 



I 1 



o 




o 




□ 


(cq Sq)p 


A 


(cq cq)^ 



FIG. 30. Eigenvector components of DD in the full basis. 



obtain the same qualitative picture for the pseudoscalar. 

For the axial vector the situation is different. We dis- 
play the mixed and unmixed effective masses in Fig. [28l 
Note that in this case the mixed D*D mass slightly in- 
creases, relative to the unmixed result. The correspond- 
ing effective overlaps are displayed in Figs. [29l[30l and[3T] 



TABLE X. Eigenvector components in the full basis. 





(cc)p 


(cc)n 


(cggc)p 


{cqqc)n 


Vc 


0.87 (5) 


-0.03 (2) 


-0.02 (1) 


-0.50 (7) 


DiD' 


0.14 (2) 


0.02 (2) 


-0.95(15) 


0.29 (4) 


J/* 


0.91 (7) 


-0.05 (2) 


0.16 (2) 


0.38(11) 


DiD 


0.14(11) 


0.07 (2) 


-0.32 (2) 


0.93 (8) 




0.41 (4) 


0.72 (3) 


-0.23 (3) 


-0.51 (4) 


DD* 


0.63 (4) 


-0.23 (3) 


-0.73 (4) 


0.12 (3) 


X'cl 


-0.55 (6) 


0.53 (5) 


-0.49 (5) 


0.41 (6) 



The fitted results on the operator overlaps are summa- 
rized in Table [X] The relative probability of creating the 
Xci by a four-quark operator is 0.29(5). For the first ex- 
citation that we identify as a D*D molecule it is 0.53(7) 
and for the x'ci 0.36(10): all these states appear to 
undergo strong mixing. 

This casts doubt onto the validity of the mixing model 
Eq. (j43p . Setting this problem aside for the moment, we 
define unmixed Xci wavefunctions, projecting onto the cc 
components of the space spanned by our operator basis, 



IXcl 



, un 



(47) 



and similarly for the excitation, x'cii while for the D*D 
state we can define the projection onto its four-quark 
components, jZ)*!?)"". Of physical interest are the over- 
laps between these idealized unmixed states and the re- 
spective physical states. These can be obtained by com- 
puting. 



D'D 
(cc)p 



D*D 
(ec)„"(ec)„ 



(48) 
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FIG. 31. Eigenvector components of x'ci in the full basis. 



The resulting probabilities read as follows, 
\{D*D\xcir? = 0.01(1) , \{xci\D*D)n^ 



Kx'cllXcirP = |(Xcl|Xel>""P= 0.01(1) 



|(D*Z?|x^i)""P= 0.22(5), \{xci\D*D)n' 
while for the normalizations we obtain, 

|(Xci|Xci)""P = 0.47(7) , \{D*D\D*Dr? 
l(XeilXci)""l' = 0.34(5). 



0.01(1), 
(49) 
0.41(7) 



0.30(6), 
(50) 



Due to cancellations the groundstate axialvector char- 
monium does not actively participate in the mixing while 
the radial excitation and the molecular state strongly mix 
with each other. However, the truncation of the mixing 
model at 0{X) is only justifiable on a qualitative level, 
as is obvious from |(Xci|Xci)""P < 1- 



VI. SUMMARY AND OUTLOOK 

We introduced the tools necessary to study the mix- 
ing of standard charmonium states with states created 
by four-quark operators and with light quark flavour- 
singlet states. Of particular importance was the use 
of the variational generalized eigenvalue method as well 
as of improved stochastic all-to-all propagator methods. 
We introduced the staggered spin partitioning (SSP) and 
the recursive noise subtraction (RNS) methods (see also 
Ref. [4^). We also made use of the hopping parame- 
ter expansion (HPE) subtraction method ^5^] and of the 
truncated solver method (TSM) [53, [111. 

Our spin-averaged charmonium spectrum agrees fairly 
well with the experimental data. However, due to our 
unphysically heavy sea quark masses with light pseu- 
doscalar masses ranging from 1 GeV down to 280 MeV, 
due to the fact that we are simulating with np = 2 sea 
quarks only and possibly due to the use of somewhat 
coarse lattice spacings of a « 0.115 fm and a « 0.077 fm. 



we significantly underestimate finestructure splittings. 
The ratio of the 25" finestructure splitting over the IS* 
splitting from which we expect some of the systematics to 
cancel comes out significantly larger than in experiment. 
One reason may be a distortion of the radial excitations 
due to their proximity to open charm thresholds which 
lie higher in our simulations than in the real world. 

The lowest spin-exotic state is a 1 vector with a 
mass of 4.15(5) GeV where the error is statistical only. 
The next highest such state can be found at 4.61(22) GeV 
with quantum numbers 2^ . We interpret the small mass 
differences that we find with respect to radial excitations 
in these sectors as evidence of a charm- ant icharm-gluon 
hybrid nature of these states. In other J^'^ sectors we do 
not see such evidence. At least for masses below 3.8 GeV 
we do not detect any mixing between S- and Z?-waves, 
indicating that L is a relatively good quantum number 
for this mass range. This conclusion is also supported by 
examining the spatial structure of the optimized creation 
operators that we employ. 

We realize that, to exclude mass shifts that are due 
to the flavour-singlet nature of charmonium states, it is 
not sufficient just to incorporate quark line disconnected 
charm annihilation diagrams but that we also have to 
consider mixing with light fiavour-singlet states. How- 
ever, within errors of less than one per mille, we do 
not detect any light quark contributions to charmonium 
wavefunctions and vice versa. Moreover, within statisti- 
cal errors of 24 MeV we do not find any significant mass 



shift: 



-.mixed 



unmixed 



11(24) MeV. Clearly, in fu- 



ture studies we will aim at reducing this error. 

We then moved on to investigate the binding between 
pairs of D and anti-il' mesons in the pseudoscalar, vec- 
tor and axialvector sectors. Only the axialvector channel 
was clearly attractive, however, we emphasize that the 
volume scaling still needs to be investigated for definite 
conclusions. Subsequently, the mixing between isoscalar 
charmonium states created by two- and four-quark op- 
erators was investigated. Within the vector and pseu- 
doscalar sector, at a light quark mass value that is four 
times as large as the physical one, these effects exist but 
they are small. However, in the axialvector channel the 
state that is dominated by the radial charmonium exci- 
tation strongly couples to the D*D molecular state and 
vice versa. This is very interesting in view of the nature 
of the experimental X{3872) state. 

We plan to apply the methods that we developed and 
tested here in high precision studies of charmonium states 
with np = 2 + 1 sea quarks of different masses on various 
volumes and lattice spacings, within a large collabora- 
tion. First results of these systematic inve stig ations were 
presented at the Lattice 2011 conference |37| . 
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